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Report is made of a systematic scaling study of the finite-temperature chiral phase transition of two-flavor QCD 
with the Kogut-Susskind quark action based on simulations on L"^ x 4 (Z/=8, 12 and 16) lattices at the quark mass 
of niq = 0.075, 0.0375, 0.02 and 0.01. Our finite-size data show that a phase transition is absent for rUq > 0.02, 
and quite likely also at rUq = 0.01. The scaling behavior of susceptibilities as a function of niq is consistent with 
a second-order transition at iriq = 0. However, the exponents deviate from the 0(2) or 0(4) values theoretically 
expected. 



1. Introduction 

Study of full QCD thermodynamics with the 
Kogut-Susskind quark action has been pursued 
over a number of years. A basic question for this 
system is the order of chiral phase transition for 
light quarks. For the case of two flavors, this 
question was examined by finite-size scaling stud- 
ies carried out around 1989-1990||l|,^. On lattices 
with the temporal size iVt = 4 and the quark mass 
in the range niq = 0.025 — 0.01, it was found that 
the peak height of susceptibilities increases up to 
a spatial lattice size L = 12, but stays constant 
within errors between L — 12 and 16. The con- 
clusion then was that a phase transition is absent 
down to rUq « 0.01, which was thought consis- 
tent with the transition being of second-order at 
niq = as suggested by the sigma model analy- 
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A more detailed study based on universality 
argument was recently attempted [||j5). Critical 
exponents were extracted from the quark mass 
dependence of the critical coupling and the peak 
height of various susceptibilities on an 8^ x 4 lat- 
tice with TOg=0.075, 0.0375 and 0.02. It was 
found that the magnetic exponent is in reasonable 
agreement with that of the 0(4) spin model ex- 
pected from universality arguments while the 
thermal exponent shows a sizable deviation from 
the 0(4) value. 

We have attempted to systematically extend 
the previous studies both regarding the spatial 
volume dependence and the quark mass depen- 
dence to further examine the universality nature 
of the transition. For this purpose we have car- 
ried out simulations on lattices of spatial size 
L = 8,12 and 16 at the quark mass of mq = 
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0.075,0.0375,0.02 and 0.01 in lattice units. In 
this article we report on results of scaling analy- 
ses based on these runsQ. Studies similar to ours 
are being carried out by other groups (7|| 

2. Simulation 

The full QCD system we study is defined by 
the partition function 

Z = Jl[dUieMSg)det{D)^f/^ (1) 

with Sg the standard single-plaquette gauge ac- 
tion, and D the Kogut-Susskind quark operator. 
Simulations are made on i"^ x 4 lattices with L — 
8, 12 and 16. For the quark mass m^, we employ 
nig = 0.075, 0.0375, 0.02 and 0.01 for each spa- 
tial lattice size L. The hybrid R algorithm]^ is 
adopted to update gauge configurations. In Ta- 
ble |l|, we list the values of (3 where our runs are 
made. To control systematic errors of the algo- 
rithm, we choose the molecular dynamics step size 
to be St « mg/2 as listed in Table |l|. For each 
run, 10000 trajectories of unit length are gener- 
ated starting from an ordered configuration. Two 
runs are made for nig = 0.01 on a 12^ x 4 lattice 
since the first run at /3 = 5.266 appears to be 
predominantly in the low-temperature phase (see 
Fig. ^ below). Critical exponents we obtain for 
L = 12 using two runs separately, however, agree 
within our statistical errors. We therefore show 
results obtained with the first run in this article. 

Inversion of the quark operator is made with 
the conjugate gradient algorithm, reducing the 
number of fioating point operations by half 
through the even-odd decimation procedure. The 
stopping condition for the even part of the source 
vector be is ^/\\be - {DWx)e\\^/3V < 10^6 ^ith 
V the space-time volume V = x A. 

Observables are calculated at every trajectory. 
For computing average values of observables we 
discard the initial 2000 trajectories of each run. 
The errors are estimated by the Jackknife method 
with a bin size of 800 trajectories. Values of ob- 
servables in the region of /3 around the simulation 
point are evaluated by the standard reweighting 
technique [p^. 

The numerical calculations have been per- 
formed on the Fujitsu VPP500/80 supercomputer 



Table 1 

Parameters of our runs. 



L 


nig 


= 0.075 


0.0375 


0.02 


0.01 




St 


= 0.05 


0.02 


0.01 


0.005 


8 


P 


= 5.35 


5.306 


5.282 


5.266 


12 




5.348 


5.306 


5.282 


5.266 












5.2665 


16 




5.345 


5.306 


5.282 


5.266 



at KEK. 



3. Observables 



In the course of our simulation, we measure the 
following susceptibilities: 



Xm 


= V 






(2) 


XtJ 






Xt.i 


- y 


\{{^i^)P,)-{Hj){p,)] 


(4) 


Xe.J 


= V 




(5) 


Xe,t 


= V[{{^Do^)P^)-{^Do^J){P,)], (6) 


Xe,ij 
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(7) 



where Dq denotes the temporal component of the 
Dirac operator, i,j = a, r, and Pa^r the spatial 
and temporal plaquette. 

Calculation of the fermionic susceptibilities 
Xm, Xt,f and Xej is non-trivial because of the 
presence of disconnected double quark loop con- 
tributions. We use the volume source method 
without gauge fixing |l^ to evaluate these suscep- 
tibilities. 

Let us illustrate our procedure for Xm- 

Per- 
forming quark contractions and correcting for the 
fiavor factor arising from the four-flavor nature of 
the Kogut-Susskind quark field, we find 



Xm 
Xdisc 



\disc 
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Xc 
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(8) 
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Let us define the quark propagator for unit source 
placed at every space-time site with a given color 
bhy 



Gl'^EiD-'C- (11) 
y 

From G"''', we calculate four quantities 0,;(i — 
1, 4) defined by 



x,y a,b 

^2 = EE^^''^^'"' 

x,y a,b 

^4 = EE^^'"^''"- 

X a,b 

It is then straightforward to show that 



(12) 
(13) 
(14) 
(15) 



(16) 



(17) 



up to terms which are gauge non- invariant, and 
hence do not contribute to the average over gauge 
configurations. We note that Oi contains con- 
nected contributions in addition to the dominant 
disconnected double quark loop contribution, and 
vice versa for 02- The terms O3 and O4 rep- 
resent contact contributions in which the source 
and sink points of quark coincide. 

4. Finite-size scaling analysis 

We start examination of our data with an anal- 
ysis of spatial volume dependence of susceptibil- 
ities for each quark mass. Let Xm""^ the peak 
height of Xm as a function of /? evaluated with the 
reweighting technique. In Fig. |l| we plot the peak 
height Xm°^ a- function of the spatial volume. 

For the heavier quark masses of niq ~ 0.075 
and 0.0375 the peak height increases little over 
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Figure 1. Peak height of the chiral susceptibility 
Xm as a function of spatial volume L^. For L = 12 
and TTig = 0.01 the upper point is from the run at 
P = 5.266 and the lower one from (3 — 5.2665. 

the sizes L = 8 — 16, clearly showing that a phase 
transition is absent for these masses. For rriq = 
0.02 an increase of the peak height is seen between 
L = 8 and 12. The increase, however does not 
continue beyond L =12; the peak height stays 
constant within errors between L = 12 and 16. 
We conclude absence of a phase transition also 
for ruq — 0.02 confirming the previous workj^,|^. 

For the lightest quark mass — 0.01 em- 
ployed in our simulation, we observe a large in- 
crease of the peak height between L — 8 and 12. 
Furthermore, the increase continues up to L = 16. 
The size dependence is consistent with a linear 
behavior in spatial volume, which one expects for 
a first-order phase transition. Other susceptibili- 
ties exhibit a similar size dependence as the quark 
mass is decreased from rUq — 0.075 to 0.01. 

This behavior contrasts with the results of 
a previous study which found that the peak 
height of susceptibilities for i = 16 stays consis- 
tent with those for L = 12 at niq w 0.01 |l). It 
is likely that a smaller statistics (2500 trajecto- 
ries||] as compared to 10000 employed here) led 
to an underestimate of susceptibilities in ref . |^ . 

An important question is whether a linear in- 
crease seen in Fig. |] could be regarded as evidence 
for a first-order phase transition at rUq = 0.01. 
We think that this is not so for several reasons. 
Looking at the time histories of the chiral order 
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Figure 2. Time history of the chiral order param- 
eter ipip for the runs with mq~O.Ol. 

parameter ipijj shown in Fig. ||, we observe an 
apparent flip-flop behavior between two different 
values of i/j^^ for L = 8. However, the time his- 
tories for L = 12 and 16 are more dominated by 
irregular patterns, and the width of fluctuation is 
smaller. These features are also reflected in the 
histograms. While we clearly see a double-peak 
distribution for i = 8, it is less evident for L — 12 
and barely visible for L = 16. Furthermore, the 
width of the distribution is smaller for larger lat- 
tice sizes and the distance between the position 
of two possible peaks is narrower. 

These observations suggest the possibility that 
the increase of the peak height seen for niq = 0.01 
up to L = 16 is a transient phenomenon due to 
insufflcient spatial volume, similar to an increase 
observed between L — 8 and 12 for rriq — 0.02. 
In order to check this point, we attempt to nor- 
malize the lattice volume by a relevant length 
scale, which we take to be the pion correlation 



length — l/'ni-n at zero temperature. Using a 
parametrization of available data for pion mass 
as a function of (3 and niq by the MILC CoUab- 
oration||lJ], we flnd ~ 3.0 for niq = 0.02 and 
4.4 for nia = 0.01. Hence the size L = 8 



0.02 roughly corresponds to L = 12 for 



for vrLq 

niq = 0.01, and i = 12 to L = 16. When com- 
pared in this correspondence the histograms for 
niq — 0.02 and 0.01 are similar in shape. It is 
quite possible that the peak height for niq = 0.01 
levels off if measured on a larger lattice, e.g., 
L = 24. 

While a definitive conclusion has to await sim- 
ulations on larger spatial sizes, we think it likely 
that a first-order phase transition is absent also 



at rrin 



0.01. 



5. Analysis of quark mass dependence 

5.1. Scaling lavifs and exponents 

We have seen in the previous section that the 
spatial volume dependence of our data do not 
show clear evidence of a phase transition down 
to niq = 0.01. In the present section we assume 
that the two-flavor chiral transition is of second- 
order occurring at niq — 0. Various scaling laws 
follow from this assumption for the quark mass 
dependence of the susceptibilities, from which we 
can extract information about critical exponents. 

For a given quark mass niq, let g^'^{niq) be the 
peak position of the chiral susceptibility Xm as 
a function of the coupling constant and let 
Xm°'^{^q) be the peak height. These quantities 
are expected to scale toward ruq = as 



^{niq) gc^{0)+Canif,<' 

max ( 



9c V-q 
Xm ("^9 



(18) 
(19) 



The peak height of other susceptibilities similarly 

scales as 

xTri'^M) - ct,, m--*-, i^f,a,T (20) 

X"r("ig) = Ce.. m---, z = /,a,r (21) 

x"rK) = ^e,^,m^'^-\ z,j=a,r (22) 

We note that Xt,i form three singular parts of the 
thermal susceptibihty xt=l^ [((V'V') e) ~ (V'V')(e)] 
with e the energy density, and Xe,i and Xe.zj 
form six singular parts of the speciflc heat C = 
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Table 2 

Critical exponents extracted by fits of critical cou- 
pling and peak height of susceptibilities for fixed 
spatial size L as compared to 0(2), 0(4)[l4|-p^ 
and mean-field (MF) values. 





0(2) 0(4) MF L = 8 


L = 12 


L = 16 




0.60 


0.54 


2/3 0.70(11) 0.74(6) 


0.64(5) 




0.79 


0.79 


2/3 0.70(4) 


0.99(8) 


1.03(9) 


Zt 


0.39 


0.33 


1/3 






Ztj 






0.42(5) 


0.75(9) 


0.78(10) 


Zt,a 






0.47(5) 


0.81(10) 0.82(12) 


Zt,T 






0.47(5) 


0.81(9) 


0.83(12) 


Ze 


-0.01 


-0.13 









ZeJ 






0.21(4) 


0.28(7) 


0.38(7) 


^e,a 






0.25(6) 


0.56(11) 


0.58(13) 


^e,T 






0.22(6) 


0.52(10) 


0.55(12) 








0.18(5) 


0.46(8) 


0.43(10) 








0.20(5) 


0.51(9) 


0.50(12) 








0.19(5) 


0.48(9) 


0.47(11) 



V [(e^) — (e)^] . The leading exponents zt and 
for xt and C are then given by zt — Max{zt_i} 
and Ze = Max{ze,i, Ze,ij]- 

For a second-order chiral phase transition, we 
expect the exponents to be expressed in terms of 
the thermal and magnetic exponents t/t and y™; 



Zg = Vt/Vh, (23) 

Zm = 2-d/yh, (24) 

Zt = l + yt/yh~d/yh, (25) 

Ze = 2yt/yh-d/yh- (26) 

Therefore two relations exist among the four ex- 
ponents Zg, Zm, Zt and Ze, which we take to be 

Zg + z„^ = Zt + 1 (27) 
2zt - z,n = Ze- (28) 



The natural values to expect for the exponents 
are those of 0{2) corresponding to exact U{1) 
symmetry of the Kogut-Susskind quark action at 
finite lattice spacing. Sufficiently close to the con- 
tinuum limit, we may also expect the 0(4) values 
as predicted by the effective sigma model analy- 
sis. The possibility of mean-field exponents ar- 
bitrarily close to the critical point has also been 
suggested p^. 
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Figure 3. Peak height of the chiral susceptibil- 
ity Xm as a function of niq for fixed spatial size 
L. Solid lines are fits to a single power (p^. 
Dashed line indicates the slope expected for 0(2) 
and 0(4) exponents which are very similar. 

5.2. Results for exponents 

Our results for the exponents are tabulated in 
Table ^. The exponent Zg that governs the scaling 
behavior of the critical coupling 5^7^ (m^) is ex- 
tracted from the fit of form (|l^) . We observe that 
Zg does not have a clear size dependence within 
our error of about 10%, and that the values are 
similar to 0(2), 0(4) or mean-field predictions, 
also listed in the Table, within one to two stan- 
dard deviations. 

Let us turn to the exponents determined from 
the peak height of susceptibilities. The values in 
Table || are extracted by fits employing a scaling 
behavior with a single power as given in (^^|-^2|) . 
In Fig. ^ we illustrate the fit for the quark mass 
dependence of the peak height of the chiral sus- 
ceptibility Xm- 

We observe in Table || that all the exponents 
z„i , Zt and Ze increase as we increase the spatial 
lattice size L. The value of Zm for the smallest 
size L —8 is not so different from the 0(2) and 
0(4) values. It deviates from the theoretical pre- 
diction for L — 12 and 16, however, and takes a 
value about 20 % larger, which amounts to a two 
standard deviation difference. 

For Zf and Zg various susceptibilities defined 
in (||)-(0) generally give consistent results. We 
observe, however, a 10 — 20% larger value of zj 
compared with the theoretical prediction already 
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Figure 4. Consistency check of exponents for a 
given spatial size L; zt + 1 against Zg + Zm, and Ze 
against 2zt — z^. Lines are predictions for 0(2) 
symmetry. Values for 0(4) are similar. 

for L — ^, and the discrepancy increases to a 
factor two difference for L = 12 and 16. The 
disagreement is more apparent for the exponent 
Ze for which values in the range Ze ~ 0.5 — 0.6 are 
obtained for larger sizes in contrast to a negative 
value for the 0(2) and 0(4) theories. 

We have noted in Sec. |5.l| that the four expo- 
nents Zg,Zm,zt and z^ should satisfy two consis- 
tency equations reflecting the fact that two rele- 
vant operators govern a second-order phase tran- 
sition. In Fig. Ewe plot the two sides of the equa- 
tions (27) and (psi) using the values of exponents 
given in Table y. For zt and z^ we take an aver- 
age over operator combinations since the values 
are mutually in agreement within the error. We 
observe that the consistency is well satisfied for 
each spatial volume even though values of individ- 
ual exponents deviate from those of 0(2), 0(4) or 
mean-field theory predictions. 

We have also attempted fits allowing for a con- 
stant term in the fitting function xT"'^ — '^oi + 
ciiTO"^'. We are not able to obtain reliable fits 
taking Zi as a free parameter, since the errors of 
fitted values are too large. Fixing the exponent 
Zi to the theoretical 0(2) or 0(4) values, we find 
that the quality of fit generally worsens compared 
with the single power fit. In particular, the fit 
tends to misses the point for the smallest quark 
mass niq = 0.01 for L = 16. 
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Figure 5. Scaling function Fm{x) normalized 
as Xmig^T'mq) ■ (mg/O.Ol)^" as a function of 
X = {i6/gl{mq) - 6/5^(0)) • (m^/O.Ol)-^' for 
L = 16 with measured values Zg — 0.6447, z„i = 
1.033,6/52(0) = 5.2353. 

We are led to conclude that the exponents show 
deviation from 0(2) or 0(4) values, at least in the 
range of quark mass nfiq = 0.075 — 0.01 explored 
in our simulation. 

5.3. Results for scaling function 

For a second-order phase transition, the singu- 
lar part of the chiral susceptibility Xm(i?^,"ig) is 
expected to scale as 



(29) 



where F„i{x) is a function of scaling variable x 
which we take to be 

X = {&/gl{mq) ~ 6/g2(0)) • m'^^. (30) 

We show in Fig. |5| estimates of the scaling 
function using data for the size L — 16. Both 
F„i{x) and x are normalized by the values for 
niq = 0.01, and the measured values are employed 
for the exponents: Zg = 0.6447, z,„ = 1.033 and 
6/5^(0) — 5.2353. Given the magnitude of statis- 
tical error which increases from 10% to 20% as 
decreases from niq = 0.075 to 0.01, we find scal- 
ing with respect to quark mass to be reasonably 
satisfied. 

We have also calculated the scaling function 
Fm{x) using the 0(4) values for the exponents 
Zg = 0.538, z„ = 0.794 and the value of 6/3^(0) 



7 



obtained with a fit of g^{mq) with the 0(4) value 
forzg. We find that scahng worsens. In particular 
the curve for the smallest quark mass nig = 0.01 
is too high in this case. 

6. Conclusions 

In this article we have reported results of our 
study of the two-flavor chiral phase transition 
with the Kogut-Susskind quark action on an Nt = 
4 lattice. Our analysis of the spatial volume 
dependence of the peak height of susceptibili- 
ties confirms the absence of a phase transition 
for niq > 0.02 as reported previously ||l|,||. At 
rriq = 0.01 the peak height exhibits an almost lin- 
ear increase over the sizes L — 8 — 16 contradict- 
ing a previous workj^. We have argued, based 
on an examination of fluctuations of observables 
and a consideration of spatial volume normalized 
by the zero-temperature pion mass, that the in- 
crease is a transient phenomenon arising from an 
insufficient spatial volume. We conclude that a 
first-order transition is likely to be absent also at 
ruq = 0.01. 

We have also found that the quark mass de- 
pendence of susceptibilities is consistent with a 
second-order transition located at niq = 0; the 
critical exponents we have obtained satisfy re- 
quired consistency conditions, and the suscep- 
tibility Xm reasonably scales in terms of vari- 
able defined with the measured exponents. How- 
ever, the values of exponents themselves deviate 
from either 0(2),0(4) or mean-field theory pre- 
dictions. 

Further work is needed to elucidate the univer- 
sality nature of the two-flavor chiral phase tran- 
sition in finite-temperature QCD. 
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